home *** CD-ROM | disk | FTP | other *** search
/ The X-Philes (2nd Revision) / The X-Philes Number 1 (1995).iso / xphiles / hp48hor2 / polyfit.src < prev    next >
Text File  |  1992-01-11  |  8KB  |  228 lines

  1. %%HP: T(3)A(D)F(.);
  2. @ POLYFIT by David F. Kurth
  3. DIR
  4. POLY
  5. @ This program takes an array of N X-Y values from the stack
  6. @ and calculates the coefficients of an N-1 polynomial that
  7. @ intersect with all N input points.  The N X-Y values are most
  8. @ easily entered using the Matrix Writer Application.
  9.  
  10. @ Directory
  11. @ Checksum # B610h
  12. @ Bytes    2267.5
  13.  
  14. \<< IF DEPTH                            @ Check for an object on stack
  15.   THEN
  16.     IF DUP TYPE 3 ==                    @ There is object, check if array
  17.     THEN
  18.       DUP 'XY' STO SIZE 1 GET 'N' STO   @ Save Array and Number of points
  19.       CLLCD                             @ Clear screen for progess updates
  20.       LijCAL                            @ Calculate L coefficients
  21.       KCAL                              @ Calculate final K coefficients
  22.     ELSE EMESS                          @ Display error message
  23.     END
  24.     {CST POLY VIEWK Fx}                 @ Restore main variable order
  25.     ORDER                               @ after new variables are created
  26.   ELSE
  27.     EMESS                               @ Display error message
  28.   END
  29. \>>
  30. CST
  31. {                                       @ Custom Menu for Polyfit program
  32.   {"POLY" POLY}                         @ Main program to find coefficients
  33.   {"VIEW" VIEWK}                        @ View polynomial coefficients
  34.   {"Fx" Fx}                             @ Given X (on stack) return f(X)
  35.   {}                                    @ Nothing
  36.   {}                                    @ Nothing
  37.   { "X-Y" XY}                           @ Recall input array to stack for editing
  38. }
  39. Cab
  40. \<< @ Generate terms to sum for C(a)(b) @
  41.   IF DUP2 ==
  42.   THEN DROP2 1                          @ If a=b, the Cab=1
  43.   ELSE
  44.     0 1 \-> a b SUM Index
  45.     \<<
  46.       { 1 }                             @ Initial list
  47.  
  48.       DO  @ For all possible terms @
  49.         WHILE @ Build list up to a-b terms @
  50.           Index a b - <                 @ Until we hit last term
  51.         REPEAT
  52.           DUP Index GET 1 + 1 \->LIST +  @ Add next term to list
  53.           Index 1 + 'Index' STO
  54.         END
  55.  
  56.         DO  @ Increment last terms until limit @
  57.           FETCH                         @ Get X values and multiply
  58.           SUM + 'SUM' STO               @ Add to sum and save
  59.           DUP Index GET 1 + Index SWAP
  60.           DUP 4 ROLLD PUT               @ Increment last term and save
  61.         UNTIL SWAP b Index + >          @ Until we go over b+Index limit
  62.         END
  63.  
  64.         IF Index 1 >  @ Increment previous terms @
  65.         THEN                            @ If possible
  66.           DO
  67.             Index 1 - 'Index' STO       @ Decrement index
  68.             1 Index SUB                 @ Shorten list
  69.             DUP Index GET 1 + Index SWAP
  70.             DUP 4 ROLLD PUT             @ Increment this term and save
  71.           UNTIL
  72.             SWAP b Index + \<=          @ This term at limit
  73.             Index 1 == OR               @ or this is last term
  74.           END
  75.         END
  76.       UNTIL
  77.         DUP Index GET b Index + >
  78.         Index 1 == AND
  79.       END
  80.       DROP                              @ Remove last list from stack
  81.       SUM
  82.       -1 a b - ^ *                      @ Fix sign of product
  83.     \>>
  84.   END                                   @ If a=b
  85. \>>
  86. EMESS
  87. \<< @ Display error message if no stack value, or value is not an array @
  88.   CLLCD
  89.   "Input array not found." 1 DISP
  90.   "Use Matrix Writer to " 3 DISP
  91.   "enter your X Y values" 4 DISP
  92.   "into an array on the" 5 DISP
  93.   "stack.  Then start" 6 DISP
  94.   "POLY again." 7 DISP
  95.   7 FREEZE
  96. \>>
  97. EQ
  98. Y.EQ                                    @ For PLOTR to plot f(x)
  99. FETCH
  100. \<< @ Use Index list on stack to fetch X values and multiply @
  101.   DUP SIZE 1 \-> lst n product          @ Save term list, size, and product
  102.   \<<
  103.     XY
  104.     1 n
  105.     FOR i
  106.       DUP lst i GET 1 XYindex GET       @ Get Xi term
  107.       product * 'product' STO           @ Multiply and save result
  108.     NEXT
  109.     DROP                                @ Drop XY array
  110.     lst product                         @ Restore term list and product
  111.   \>>
  112. \>>
  113. Fx
  114. \<<  @ Calculate F(x) using K coefficients and x value from stack @
  115.   K N GET                               @ First coefficient
  116.   N 1 - 1
  117.   FOR i                                 @ For i = N-1 to 1
  118.     OVER *                              @ x times
  119.     K i GET +                           @ Add in next coefficient
  120.     -1
  121.   STEP
  122.   SWAP DROP                             @ Remove x value
  123. \>>
  124. KCAL
  125. \<< @ Calculate K coefficients (f(x) = Kn*x^n-1 + ... K2*x + K1) @
  126.   1 N
  127.   START                                 @ Create K array with N zeroes
  128.     0
  129.   NEXT
  130.   N \->ARRY 'K' STO                     @ Save K array
  131.  
  132.   1 N
  133.   FOR i                                 @ i=1 to N
  134.     "DOING K" STD i \->STR + 1 DISP     @ Display Ki value
  135.     0                                   @ Zero sum
  136.     i 1 - N 1 -
  137.     FOR j                               @ j=i-1 to N-1
  138.       @ Ki = Ki + Lj1 * C(j)(i-1) @
  139.       L j 1 Lindex GET                  @ Get L term from array
  140.       j i 1 - Cab                       @ Get Cab coefficient
  141.       * +                               @ Sum to Ki
  142.     NEXT
  143.     K i 3 ROLL PUT 'K' STO              @ Add K value to array
  144.     "K" i \->STR + "=" + K i GET \->STR + 3 DISP
  145.   NEXT
  146. \>>
  147. LijCAL
  148. @ Calculate triangular table of L values @
  149. \<<
  150.   "DOING L:" 1 DISP
  151.   1 TMAX                                @ Create L values array
  152.   START                                 @ With TMAX zeroes
  153.     0
  154.   NEXT
  155.   TMAX \->ARRY                          @ L array initialized
  156.   1 N                                   @ Set Lj,1 = Yj
  157.   FOR j                                 @ For i=0, j = 1 to N
  158.                                         @ L array already on stack
  159.     j                                   @ L index value to place Y value
  160.     XY j 2 * GET                        @ Y value for L array
  161.     PUT                                 @ Y value now in L array
  162.   NEXT
  163.                                         @ Initialized L array left on stack
  164.   1 N 1 -
  165.   FOR i                                 @ For i = 1 to N-1
  166.     1 N i -
  167.     FOR j                               @ For j = 1 to N-i
  168.       DUP DUP i 1 - j 1 + Lindex GET    @ Get L i-1,j+1 value
  169.       SWAP i 1 - j Lindex GET -         @ Get L i-1,j value and subtract
  170.       XY DUP i j + 1 XYindex GET        @ Get X i+j value
  171.       SWAP j 1 XYindex GET -            @ Get X j value and subtract
  172.       / i j Lindex SWAP PUT             @ Save new L i,j into array
  173.     NEXT
  174.   NEXT
  175.   'L' STO                               @ Save final array
  176. \>>
  177. Lindex
  178. @ Return index into L array of i,j element @
  179. \<<
  180.   \-> i j                               @ Save index values
  181.   \<<
  182.     i 1 - N i 2 / - * N + j +           @ index=(i-1)(N-i/2)+j+N
  183.   \>>
  184. \>>
  185. PPAR
  186. {
  187.   (-6.5,-3.1)                           @ X range
  188.   (6.5,3.2)                             @ Y range
  189.   X                                     @ Independent Variable
  190.   0                                     @ Resolution (pixel in each column)
  191.   (0,0)                                 @ Axes intersection
  192.   FUNCTION                              @ Function type
  193.   Y                                     @ Dependent Variable
  194. }
  195. TMAX
  196. @ Maximum Number of elements in the L triangular array @
  197. \<<
  198.   N DUP 1 + * 2 /                       @ TMAX = N(N+1)/2
  199. \>>
  200. VIEWK
  201. \<< @ Display the K coefficients of the resulting polynomial expression @
  202.   STD CLLCD "Hit VIEW to cont..." 1 DISP
  203.   N 1
  204.   FOR i                                 @ For each coefficient
  205.     "K" i \->STR + "=" +
  206.     K i GET \->STR + "*x^" + i 1 -
  207.     \->STR + 3 DISP 0 WAIT DROP
  208.     -1
  209.   STEP
  210.   CLLCD                                 @ Clear display to let user know
  211.                                         @ the program is still working
  212.   KILL                                  @ Get rid of HALT annunciator
  213. \>>
  214. XYindex
  215. @ Return index into XY array of i,j element @
  216. \<<
  217.   \-> i j                               @ Save index values
  218.   \<<
  219.     i 1 - 2 * j +                       @ index=2(i-1)+j
  220.   \>>
  221. \>>
  222. Y.EQ
  223. \<< @ Fx function requires X value from stack, not compatible with PLOTR. @
  224.   @ This program takes 'X' from PLOTR function, calls Fx, and @
  225.   @ returns with f(x) value to satisfy PLOTR requirements. @
  226.   X Fx
  227. \>>
  228.